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Forward 

The proposed research is to develop Gaussian random fields methods to study fork-join net¬ 
works (FJNs) with synchronization constraints. FJNs arise from many military operations, e.g., 
Army force deployment and counter-terrorism, where commands come from one or multiple types 
of operations and each operation requires multiple parallel and/or sequential tasks to be processed 
in service stations with multiple servers, and to be rejoined for further processing with synchro¬ 
nization constraints, e.g., non-exchangeability. In this research, we focus on the non-exchangeable 
synchronization constraint, which requires that tasks can only be synchronized only if all tasks of 
the same job are completed. The main mathematical challenge lies in the resequencing of arrival 
orders after service completion at each station, which requires an infinite dimensional state space 
to track the status of all parallel tasks for each job. That was an extremely difficult open problem. 

We have developed a novel method using multiparameter sequential empirical processes driven 
by service vectors of parallel tasks of each job to describe the system dynamics of FJNs. This 
research has produced two research papers, focusing on a single class FJN in two asymptotic regimes, 
where the arrival rate of jobs and the number of servers in each station get large appropriately. We 
consider the number of tasks in each waiting buffer for synchronization, jointly with the number 
of tasks in each parallel service station and the number of synchronized jobs. In the first paper, 
we consider the quality-driven regime, and show that all the limiting processes are functionals of 
two independent processes - the limiting arrival process and a generalized Kiefer process driven 
by the service vector of each job. We characterize the transient and stationary distributions of 
the limiting processes. In the second paper, we consider the quality-and-efficiency-driven regime 
(Halfin-Whitt regime), and show that all the limit processes in the functional central limit theorem 
are also characterized via functionals of the initial limit quantities, the arrival limit process and 
a generalized multiparameter Kiefer process driven by the service vectors. This new framework 
is being further generalized to analyze fork-join networks with multiple classes of jobs, and study 
control, reliability and provisioning problems. 
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1 Statement of the Research Problem 


Fork-join networks consist of a set of service stations that serve job requests simultaneously and 
sequentially according to pre-designated deterministic precedence constraints. Such networks have 
many applications in manufacturing and telecommunications [4, 16, 25, 26, 27, 43, 53, 36, 37, 49], 
patient flow analysis in healthcare [22, 1, 2, 57, 58], parallel computing [47, 52, 51, 32], military 
deployment operations [24, 56], and law enforcement systems [29]. Two types of synchronization 
constraints are of particular interest. One is called exchangeable synchronization (ES) in which 
tasks are not tagged with a particular job and can be synchronized for a service completion once 
the necessary tasks are completed. This type of synchronization constraint is often used in manufac¬ 
turing systems; for example, in many assembly systems, different parts of a product are processed 
at separate workstations or plant locations and a product will be assembled once all of its neces¬ 
sary parts are completed. In this case, the parts are not tagged with a particular product, since 
they are standardized for the same type of product. The second type is called non-exchangeable 
synchronization (NES). Tasks are tagged with a particular job and can only be synchronized when 
all the parallel tasks of the same job are completed. 


Fort; Service Stations Unsynchronized Queues 



Figure 1: A fundamental fork-join network 

Fork-join networks with NES are used in many applications, including healthcare systems, 
parallel computing, MapReducing scheduling (e.g., large-scale parallel Web search), disassembly 
and reassembly systems in manufacturing and so on. In patient flows of hospitals [1, 2, 22, 57, 58], 
the treatment and discharge processes are typical examples of fork-join networks with NES: a 
patient must have all test results ready before a doctor examination and these tests are conducted 
in different units/laboratories and can never be mixed; a patient, after the discharge decision is 
made, must wait for necessary procedures, pharmacy, transportation, etc., before being physically 
discharged. In MapReduce scheduling [11, 32, 51, 54], jobs are processed in two phases: in the 
map phase, a large-scale data input (e.g., Web processing data) is distributed into individual 
computation nodes, and each node processes one block of input data, and after the execution of all 
blocks of the same data input, they will be joined as an output in the reduce phase. 

Despite the vast appealing applications of such networks, very little has been known about their 
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behaviors in the many-server heavy-traffic regimes. We start considering a fundamental fork-join 
network model with a single class of jobs and NES, where each arriving job is forked into several 
parallel tasks upon arrival and each of the tasks is processed in parallel at a dedicated service 
station with multiple servers under the non-idling FCFS discipline. Upon service completion, each 
task will join a buffer associated with its service station, and wait for synchronization, such that 
each job is synchronized only if all of its tasks have been completed. Figure 1 depicts such a 
model. In this model, in addition to the service dynamics, we are interested in the waiting buffer 
dynamics for synchronization. One important performance measure is the response time of a job, 
namely, the time from arrival to synchronization. The response time may also include the time 
required for the synchronization process, but we do not consider that in this work. Thus, the 
response time includes two delays, waiting time for service and waiting time for synchronization. 
Since each service station can be regarded as a separate many-server queue, the waiting time for 
service has been well understood. However, the waiting time for synchronization, which is our focus 
in this paper, has not been studied. Specifically, we investigate the waiting buffer dynamics for 
synchronization jointly with the service dynamics. 

The main mathematical challenge lies in the resequencing of the arrival orders after service 
completion at each service station, due to the randomness of the service times and the multi-server 
setting. When there is a single server in each of the parallel service station and the service discipline 
is FCFS, the service completion order is preserved to be the same as the arrival order of tasks in 
each service station, so that the two types of synchronization constraints are equivalent. However, 
the arrival order of tasks in each service station can be resequenced at the service completion 
epochs when the number of servers in a service station is larger than one or the service discipline 
is not FCFS. Resequencing has been one of the most difficult obstacles in the study of fork-join 
networks. Some limited work has been dedicated to the study of such challenging problems. For 
example, substantial efforts were dedicated to the study of the max-plus recursions [21, 3, 12]. 
More recently, Atar et al. [2] have studied a fork-join network with single-server service stations 
where tasks may reenter for service at some service stations in a Bernoulli mechanism so that 
the arrival orders of tasks at each service station are resequenced after service completion. They 
show that under a priority discipline, the system dynamics with NES is asymptotically equivalent 
to that with ES in the conventional (single-server) heavy-traffic regime. For a Markovian fork- 
join network with multiple servers, Zviran [58] shows that the system dynamics with NES is also 
asymptotically equivalent to that with ES in the conventional heavy-traffic regime. However, the 
two types of synchronization constraints lead to very different system dynamics when the service 
stations have many parallel servers in the Halfin-Whitt regime, as conjectured in [2, 58]. To the 
best of our knowledge, our work is the first to tackle the resequencing problem in non-Markovian 
fork-join networks with NES and multiple-server service stations in the many-server heavy-traffic 
regimes. We will consider both cases when each service station is operating in the quality-driven 
(QD) regime, or in the quality-and-efficiency-driven (QED, Halfin-Whitt) regimes. 

When all the service stations operate in the QD regime, this is equivalent to a model which has 
infinite numbers of servers at all service stations asymptotically. To describe the system dynamics, 
we can start with a graphical representation as shown in Figure 2(a) for a system of two parallel 
tasks. At each job’s arrival epoch, we mark the arrival time on the horizontal line (rc-axis) and 
the service times of all parallel tasks on the vertical line (y-axis). At each time t , by drawing 
a negative forty-five degree line, we can count the numbers of tasks in each service station and 
each waiting buffer for synchronization. When the arrival process is Poisson, we can apply Poisson 
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(a) The QD Regime 



Figure 2: Graphical representations of the system dynamics in the QD and QED regimes 


random measure theory, similarly as in the “physics” of M/GI/oo queues [14]. It can be shown 
that at each time t, the numbers of tasks in each service station and each waiting buffer for 
synchronization all have Poisson distributions and their parameter values and covariances can also 
be obtained; see Proposition 2.1. However, when the arrival process is more general, this Poisson 
random measure approach does not work, and we cannot obtain the exact distributions for these 
performance measures. Thus, we consider heavy-traffic approximations of the system dynamics 
when the arrival rate is relatively large. For that, the graphical representation in Figure 2(a) also 
plays an important role; see the system’s dynamic equations in §2. 

Here we develop a new approach to describe the system dynamics. Both the service dynamics 
and the waiting buffer dynamics for synchronization are represented as functionals of the mul¬ 
tiparameter sequential empirical process driven by the service vector of all parallel tasks. Their 
diffusion-scaled processes converge weakly to limit processes that can be all represented as function¬ 
als of two independent processes - the limiting arrival process and the multiparameter generalized 
Kiefer process driven by the service vector. When the limiting arrival process is Brownian motion, 
we show that the aforementioned limiting processes are a multidimensional continuous Gaussian 
process, and thus characterize the joint transient and stationary distributions of these processes. 
We also study the impact of the correlation among the service vector upon these distributions. 

There are several advantages with this new approach. It gives a clean and elegant representation 
of the limiting processes, involving only two independent stochastic processes arising from the arrival 
and service processes. Moreover, the characterization of the limiting processes as Gaussian and their 
transient and stationary distributional properties can be easily obtained. Furthermore, this new 
approach paves the way to study the fork-join network with all the service stations operating in the 
QED regime. We believe that this new approach launches a new framework to study more general 
fork-join networks, for example, multiclass models, and when the service vectors for parallel tasks 
form a stationary and weakly dependent sequence. 

When all the service stations are in the QED regime, we exploit the delicate relationship between 
finite-server models and its corresponding infinite-server models. This was exploited to prove 
an FCLT for the GI/GI/n queue by Reed [45]. We make an important observation that the 
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multidimensional processes of the waiting buffer dynamics for synchronization and the service 
dynamics in the fork-join network can be represented through the corresponding processes in the 
infinite-server case. Thus, our results from the QD regime can be extended to establish the FCLT 
for the fork-join network in the QED regime. To illustrate, we can also use a similar graphical 
representation as in Figure 2(a) to describe the system dynamics. In particular, as shown in Figure 
2(b), we mark the entering service times of all parallel tasks for each job on the horizontal line {x- 
axis), and the service times of them on the vertical line (y-axis). However, unlike the infinite-server 
case, tasks of the same job may not enter service simultaneously. Fortunately, it is well known 
that the delay for service in the QED regime is 0{\/y/n)\ see, e.g., [45, 50]. This asymptotically 
negligible difference among entering service times helps us to establish the FCLT for the fork-join 
network in the QED regime. 

An important implication of our results is that the size of the waiting buffer for synchronization 
is of the same order as that of the total number of tasks at each service station, and thus, the waiting 
time for synchronization is of the same order as the service time, 0(1). Namely, the response time 
in the QED regime includes the delay for service 0(1/ \/n), the service time 0(1) and the delay 
for synchronization 0(1). It remains to establish the FCLT for the (virtual) waiting time process 
for synchronization. More importantly, it remains to find an optimal scheduling policy that will 
minimize the delay for synchronization in the single-class case. We believe that our methods and 
results will provide useful insights towards that direction. 

In the development of approximations to the fork-join system, we make a fundamental contri¬ 
bution to the study of multiparameter sequential empirical processes driven by random vectors. 
Sequential empirical processes driven by a sequence of random vectors (allowing for correlation 
among random variables in the vector) and their limits as generalized Kiefer processes have been 
studied in the statistics literature; see e.g., [42, 6, 8, 9, 13], but the convergence is proved in the 
space .D([0,T] fc ,R) of real-valued cadlag functions defined on [0,T] fc , k > 2, endowed with the 
generalized Skorohod Ji topology in [35] and [48]. In our setting, it is necessary to prove the con¬ 
vergence in the space _D([0, T], D([0, T] fc , R)) of function-valued cadlag functions defined on [0, T], 
endowed with the standard Skorohod J\ topology for D([0, T] k , Revalued cadlag functions. 

Literature review. Most of the literature on fork-join networks is on models with single-server 
service stations. We only give a brief summary here on relevant work in heavy traffic. These 
studies are in the conventional (single-server) heavy-traffic regime. In Varma’s dissertation [53], 
the diffusion-scaled workload processes and unsynchronized queueing processes in some fork-join 
network models with ES are shown to converge weakly to certain multi-dimensional reflected Brow¬ 
nian motions. The stationary distributions of the system response time and the processes counting 
the number of tasks in unsynchronized queues are specified by some partial differential equations 
(PDEs). Nguyen [36] shows the diffusion-scaled processes counting the queue lengths at each service 
station of a single-class fork-join network model with ES converge to a reflected Brownian motion 
in a polyhedral cone of the nonnegative orthant. Nguyen [37] discusses the difficult challenges with 
multiclass fork-join models with ES. As we have noted above, for a fork-join network with feedback 
and NES, Atar et al. [2] show that a dynamic priority discipline achieves throughput optimal¬ 
ity asymptotically in the conventional heavy-traffic regime, as a consequence of the asymptotic 
equivalence between NES and ES constraints. 

Very little work has been done for fork-join networks with multi-server service stations. Ko 
and Serfozo [25] consider a fork-join network model with a single class of Poisson arrivals and K 
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parallel service stations with multiple servers at each station and exponential service times, and 
obtain an approximation for the distribution of the system response time in equilibrium under 
the NES constraint. Dai [10] provides an exact simulation algorithm to approximate the system 
response time in equilibrium for the same Markovian model in [25] by using a “coupling from the 
past” method. Zviran [58] studies optimal control of multi-server feedforward fork-join networks 
with exponential service times in the conventional heavy-traffic regime and shows that FCFS is 
asymptotically optimal and the resequencing disruption becomes asymptotically negligible. Zaied 
[57] calculates mean offered-load functions of fork-join networks with NES and multiple processing 
stages when the arrival process is time-inhomogeneous Poisson and service times for parallel tasks 
are independent, and studies staffing of time-varying emergency departments and synchronization 
delays under Markovian assumptions. Both dissertations of Zviran [58] and Zaied [57] are motivated 
from applications in patient flow analysis. Gurvich and Ward [17] study optimal matching policies 
for a pure join model (Markovian) with multiple classes of jobs under certain matching constraints. 

This work contributes to the recent development for non-Markovian many-server queueing mod¬ 
els. We only mention those that are most relevant to our work due to the large volume of papers 
on many-server models. Krichagina and Puhalskii [28] first observe that the system dynamics of an 
infinite-server queueing model can be represented by an integral functional of a sequential empirical 
process driven by service times. They show that the diffusion-scaled processes counting the num¬ 
ber of jobs in the system can be approximated by a functional of a standard Kiefer process driven 
by service times. Pang and Whitt [39, 41] generalize that approach to establish two-paranreter 
process limits for G/G/oo queues when the service times are i.i.d. and weakly dependent, respec¬ 
tively. Reed [45] and Puhalskii and Reed [44] have observed a relationship between finite-server 
and infinite-server queues and generalized the approach in [28] to obtain the diffusion limits for 
G/GI/N queues in the Halfin-Whitt regime. Mandelbaum and Momcilovic [33] generalize the 
approach by Reed [45] to study G/GI/N + GI queues with abandonment. All these papers use 
sequential empirical processes driven by a sequence of univariate random variables. Our approach 
to study fork-join networks with NES uses multiparameter sequential empirical processes driven by 
a sequence of i.i.d. random vectors and properties of multiparameter processes and martingales. 

Notation Throughout the paper, the following notation will be used. M and (R rf and Ml, 
respectively) denote sets of real and real non-negative numbers (d-dimensional vectors, respectively, 
d > 2). Z+ is the set of non-negative integers. N denotes the set of natural numbers. For a, b G M, 
we denote a A b := min(a, b) and a V b := max(a, b). For x G M, let x + := max{x, 0} and 
x~ := — min{x,0}. For any x G M+, is used to denote the largest integer less than or equal to 
x. We use bold letter to denote a vector, e.g., x := (x\,... ,xn ) £ M^. 0 denotes the vector whose 
components are all 0. For x,y G M w , we denote x < y, x > y and x > y in the componentwise 
sense, and let x A y = (x\ A yi, ...,xn A yjsr). We use 1(A) to denote the indicator function of a 
set A. The abbreviation a.s. means almost surely. For any univariate distribution function F(-), 
we denote F c (-) = 1 — F(-). For a G and a G M + , we call A a (5) (resp . A a (5)) is a 5-grid of 
[0,ai] x [ 0 , 02 ] (resp. [0, a]), if A a (6) (resp. A a (S)) is a finite partition of [0, a{\ x [ 0 , 0 : 2 ] (resp. 
[0, a]), where each element of the partition is the rectangle [si,ti) x [S 2 T 2 ) (resp. [s,t)), satisfying 
0<s k < t k < ah for k = 1, 2 (resp. 0 < s < t), and min k=i, 2 (tk ~ s k ) > 5 (resp. t — s > 5). For 
two real-valued functions / and g, we write f(x) = 0(g(x )) if limsup^^ \ f(x)/g(x)\ < 00 . 

All random variables and processes are defined on a common probability space (D, F, P). For 
any two complete separable metric spaces <Si and S 2 , we denote <Si x S 2 as their product space, 
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endowed with the maximum metric, i.e., the maximum of two metrics on 5i and 52• 5 fc is used 
to represent fc-fold product space of any complete and separable metric space 5 for k G N. For a 
complete separable metric space 5, D([0, oo),5) denotes the space of all 5-valued cadlag functions 
on [0, 00 ), and is endowed with the Skorohod J\ topology (see, e.g., [5, 15, 55]). Denote D = 
D([0, 00 ), M). The space D([0, 00 ),B), denoted as Du, is endowed with the Skorohod J\ topology, 
that is, both inside and outside D spaces are endowed with the Skorohod J\ topology. For a 
complete separable metric space 5, the space D([0, oo) 2 ,5) is the space of all 5-valued “continuous 
from above with limits from below” functions on [0,oo) 2 , and is endowed with the same metric as 
defined by [18]. D 2 = O([0, l] 2 , M) is denoted as the space of all “continuous from above with limits 
from below” functions on the unit square [0, l] 2 in the sense of Neuhaus [35], and is endowed with 
the same metric dj$ 2 as in [35]. Weak convergence of probability measures //„ to // will be denoted 

d f 

as n n =$■ /I. For a sequence of processes {X n : n > 1} and a process X, we use notation X n A X 
to denote the convergence in finite-dimensional distributions of X n to X. 


2 The Infinite-Server Fork-Join Network Model 

2.1 Model and Assumptions 

In this section, we present a detailed description of our infinite-server fork-join network model and 
the assumptions. As shown in Figure 1, there is a single class of jobs, and each job is forked into K 
parallel tasks, K > 2. Each task is processed in a service station with multiple servers under the 
FCFS discipline. There is an infinite number of servers at each station. After service completion, 
each task will join a waiting buffer for synchronization associated with each service station, and 
when all tasks of the same job are completed, they will be synchronized and leave the system. Here 
we assume that the synchronization process takes zero amount of time. 

Let A := {A(t) : t > 0} be the arrival process of jobs with r* representing the arrival time 

of the ?' th job, i £ N. Let { 77 * : i > 1} denote the i.i.d. service time vectors of the parallel 

tasks. The joint distribution of the service time vector for the 7 th job 77 ' is F(x) := F(x 1 ,... ,xk ) 
for x k > 0, k = 1 Their marginal distributions are Ff : (x), for x > 0, k = The 

joint distribution of any two service times //* and rf k is Fj, k (xj,x k ) := P(r/) < x 3 , rf k < Xk) for 
Xj,Xk > 0, j,k = 1,..., K. Note Fjk(-,-) = Pfc(-) when j = k for j,k = 1 We denote 

F^ k (xj,x k ) := P(rfj > Xj,rf k > x k ) = 1-Fj(xj)- F k (x k ) + F jtk (xj, x k ) for xj,x k > 0 ,j,k = 

Note Fj k (-, •) = F k (-) when j = k for j,k = 1, Let 'rf m := max{r/j, ...,rf K } be the maximum 

of the components in the service vector 77 *, and F m (x) := P(rjm — x ) = F(x,...,x) for x > 0. 
(Throughout the paper, we use subscript “m” to index quantities and processes associated with 
the maximum.) The service process is assumed to be independent of the arrivals. We exclude the 
case of perfectly positively correlated parallel services since that will lead to empty waiting buffers 
for synchronization. 

Let X k := {X k (t) : t > 0} be the process counting the number of tasks in service at the service 

station k. and Y k = {Y k (t) : t > 0} be the process counting the number of tasks in the waiting 

buffer for synchronization (unsynchronized queue) after service completion at service station k, 
k = 1,..., K. Let 5 := {5(f) : t > 0} be the process counting the number of synchronized jobs and 
D k := {D k (t) : t > 0} be the process counting the number of tasks that have completed service at 



station k, k = 1 Denote X := (X\,...,Xk), Y := (Yi, Yk) and D := (D\, ..., Dk). We 

assume that the system starts empty. 

Assuming that the arrival process Aft) is Poisson with rate A, by Poisson random measure 
theory, we can easily obtain the following properties on the processes X(t), Yft) and Sft) at each 
time t. 

Proposition 2.1. If the arrival process Aft) is Poisson with rate X, then at each time t > 0, for k = 
1 Xk{t) has a Poisson distribution with rate X Ff(s)ds, Y k {t) has a Poisson distribution 

with rate A Jq(F^(s) — Ff(s))ds, and S(t) has a Poisson distribution with rate A F m (s)ds. For 
each time t > 0 and j, k = 1,..., K, 

Cov(Xj(t),X k (t )) = X f* Fj k (s,s)ds, (2.1) 

CovfYjft), Y k (t)) = A /Vi, k(s, s) - F m (s))ds, (2.2) 

Jo 

Cov(Xj(t),Y k (t )) = A [ {F k (s)-F j>k (s,s))ds. (2.3) 

Jo 

For each time t > 0 and k = 1,..., K, S(t ) is independent of X k (t) and Y k ft). When K = 2, Y\ (t) 
and Yiit) are independent for each t > 0. 

When the arrival process Aft) is general, we will obtain heavy-traffic limits for the fluid and 
diffusion scaled processes of (X, Y, S ) jointly. We will let the arrival rate grow large for the system 
to be in heavy traffic. For that, we consider a sequence of such systems indexed by n and use 
superscript n for the processes A,X ,Y, D, S, and the arrival times {r; : i > 1}, but we let the 
service times {rf : i > 1} and their distribution functions be independent of n. We make the 
following assumption on the arrival process A n . 

Assumption 1: FCLT for arrivals. There exist: (i) a continuous nondecreasing deterministic 
real-valued function a on [0, oo) with a(0) = 0 and (ii) a stochastic process A with continuous sample 
paths, such that 

A n := n ~2 ( A n — na) A in D as n —»• oo. (2-4) 

■ 

It follows from (3.6) that we have the associated FWLLN 

A n 

A n := — a in D as n —> oo. (2-5) 

n 

When the arrival process is renewal, the limit in (2.5) is aft) = At, for t > 0 and some positive 
constant A, and the limit in (3.6) is A = Xc^B a , where cf is the squared coefficient of variation 
(SCV) of an interarrival time, and B a is a standard Brownian motion (BM). 

We also make a regularity assumption on the joint service-time distribution function F(x). 

Assumption 2: Service time distributions. The joint distribution function F(x) of the 
service time vectors {q l : i £ N} is continuous. ■ 

From the graphical representation of the system dynamics in Figure 2(a), we can write, for each 
t, > 0 and k = 1,..., K , 

A n {t) 

*?(*)=£ 1 (t“ + 4 >*), (2.6) 

i= 1 
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(2.7) 


A n (t) 

Y£(t) = ^ l(r” + rf k <t and r” + rf k , > t for some k! / k) 

i =1 
A"it) 

= (W + rfk -V- 1 (tP + < < 0) 

i=l 
A-it) 

= Y1 W 7 ? + ?? m > *) “ 1 (^ + Vk > 0) , 

j=l 

A"(t) A"(t) 

S"(t) = V 1(7? + i < t) = £ !W + 4 < t, vi), (2.8) 

7=1 7=1 

A™(t) 

D n k (t)= £ 1 (r? + 4<t). (2.9) 

i=l 

The following balanced equations hold for each t > 0 and k = 1,..., K, 

Dk(t) = A n (t) — X k (t), (2.10) 

Y k n (t) = Dl(t) -S n (t). (2.11) 

As we have remarked in the introduction, by previous work on G/GI/oo queues [28], each 
individual process X k and D k (resp. S n ) can be represented by an integral of a sequential empirical 
process driven by a sequence of i.i.d. random variables {rf k : i > 1 } (resp. {rf m : i > 1 }) for each 
k = 1, Thus, Gaussian limits for the diffusion-scaled processes X k , DV and S n in heavy traffic 

for each k can be established, and as a consequence, a Gaussian limit for the diffusion-scaled process 
Y k can be obtained from those of DV and S n , k = 1,..., K. However, that approach does not give a 
characterization of the joint Gaussian distribution of the limiting processes of the diffusion-scaled 
processes (X n ,Y n ,S n ). 

We will represent all the processes X n , Y' 1 , S n as integrals of a multiparameter sequential empir¬ 
ical process K n := {K n (t,x) : t > 0,x G M+} driven by the sequence of service vectors {rf :*>!}: 


\nt\ 


K n (t,x ) := - V l(rf <x), t> 0, x € 

n ^ 


That is, we write, for t > 0 and k = 1, 


and 


X k (t) = n / / 1 (s-h x k > t)dK n (A n (s),x) 

Jo Jr* 


Y k (t) =n / (l(s + Xk < t) — l(s + Xj < t, Vj)) dK n (A n (s),x) 

Jo J Kf 


S n (t ) =n / l(s + Xj < t, Mj)dK n (A n (s),a:) 

Jo Jr* 


( 2 . 12 ) 

(2.13) 

(2.14) 

(2.15) 


The integrals in (2.13), (2.14) and (2.15) are well-defined as a Stieltjes integral for functions of 
bounded variation as integrators. 
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2.2 An FCLT for Multiparameter Sequential Empirical Processes 

We present an FCLT for multiparameter sequential empirical processes U n := {U n (t,x) : t > 0, x E 
[0,1]*} driven by a sequence of i.i.d. random vectors with uniform marginals: 

[nt.] 

U n (t,x) := —/= ^ (l(£* < x) - H(x)) , t> 0, X6[0,1] A \ (2.16) 

where for each iGN, £' := (£j,.... f l K ) is a vector of nonnegative random variables with continuous 
joint distribution function H(- ) and uniform marginals over [0,1]. 

The convergence for the processes U n (t,x ) is established in the space D([0, oo), D([0, 1 ] K , M)). 
We remark that this theorem is in the same spirit as Lemma 3.1 in [28], where an FCLT is proved 
for the two-parameter process U n (t,x ) in the univariate case in the space D([0, oo), D([0,1], R)). 
We generalize that result to the multivariate setting. 

Theorem 2.1. The multiparameter sequential empirical processes U n (t.,x ) defined in (2.16) con¬ 
verge weakly to a continuous Gaussian limit, 

U n {t,x) => U(t,x) in D([0, oo),D([0,1] A , R)) as n —>■ oo, (2-17) 

where U(t,x ) is a continuous Gaussian random field with mean function E[U(t,x)\ = 0 and covari¬ 
ance function 

Cov(U(t,x), U(s,y)) = (t A s)(H(x Ay) — H(x)H(y)), t,s> 0, x,y g [0, 1] K . 

To show the FCLT for the processes (X n ,Y n , S n ), we define the diffusion-scaled multiparameter 
sequential empirical processes K n := {K n (t,x) : t > 0,16 R;^} by 

[nt\ 

K n (t,x) := —= (l(? 7 * < x) — Fix)) , t > 0, i6lf. (2-18) 

1=1 

Theorem 2.1 can be applied to show an FCLT for the processes K n (t,,x). Define F : R K —»• [0,1] A 
with F(x) = (Fi(xi),..., Fk{xk))- By Sklar’s theorem [46], for any multivariate distribution func¬ 
tion F. there exists a unique multivariate distribution function H (called “copula”) with uni¬ 
form marginals on [0,1] such that F(x) = H(F(x)) when the marginal distribution functions Ff ., 
k = 1 are continuous. Then, K n (-,-) can be represented as a composition of U n (-,-) with 

F(-) in the second component, i.e., 

K n (t,x) = U n {t,F(x)), t> 0, xgR^. 

Thus, it follows from Theorem 2.1 that the processes K n (t,x) converge in distribution: 

K n (t,x) = U n (t,F(x)) => K(t,x) := U(t,F{x)) in D([0, oo), D#) as n —>• oo, (2.19) 
which implies that 

K n (t,x) =>■ k(t,x) := tF(x) in D([0, oo), D#) as n —> oo. (2.20) 
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2.3 FWLLN and FCLT 

We define fluid-scaled processes X n , Y n and S n by 

X n ■= l X n , Y n := -Y n , S n := -S n . (2.21) 

n n n 

The FWLLN for (X n ,Y n ,S n ) is stated in the following theorem. 

Theorem 2.2 (FWLLN). Under Assumptions 1 and 2, the fluid-scaled processes converge to 
deterministic fluid functions, 


(A n ,X n ,Y n ,S n ) => (a,X,Y,S) (2.22) 

in B 2K+2 as n —>• oo, where the limits are all deterministic functions: a is the limit in (2.5), for 
each t > 0, 


X(t) := {Xflt),...,X K {f)), X k (t) := / Ff(t - s)dd(s), for k = 1, K, (2.23) 

Jo 

Y(t) :=(Fi(f),...,y^(t)), Y k (t) := - s) - Ff(t - s))dd(s), for k = l,...,K, (2.24) 


S(t) := / F m (t — s)da(s). 

Jo 

When a(t) = A t for a constant arrival rate A > 0 and E[i)j,] < oo /or k = 1,..., A/ 
A fc (oo) := lim X k (t) = A k = 1,..., K, 

t—> OO 

Yfc(oo) := lim T fc (t) = A(A[r/^] - A^]), k = 1,..., K, 

t —/OO 

lim M = A . 

>-oo t 

We define the diffusion scaling of X n , F ra and S n by 

X n := V^(^ n -^), y n := ^n(Y n ~Y), S n := y/n(S n - S). 


(2.25) 

(2.26) 

(2.27) 

(2.28) 

(2.29) 


We will show the following FCLT for these diffusion-scaled processes. 

Theorem 2.3 (FCLT). Under Assumptions 1 and 2, the diffusion-scaled processes converge in 
distribution, 

(. A n ,K n ,X n ,Y n ,S n ) => (A, K,X,Y, S) (2.30) 

in D x D([0, oo),D a') x D 2A+1 as n —> oo, where A is the limit in (3.6), K is the limit in (2.19), 
which is independent of A, and for t > 0 and k = 1,..., K, 


X(t) 


M\(t) + M 2 (t), Miff) := (M M (f), ..., i = 1,2, 

fm - s)dA(s) = A(t) - f A{s)dF c k {t - s), 
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(2.31) 

(2.32) 



Mk,2(t) 
S(t) 
Vi (t) 

V 2 (t) 

Y(t) 

Zk,i (t) 
%k,2 (t) 


l(s + Xk > t)dK(a(s),x) = — / / 1 (s + x k < t)dK(a(s),x), (2.33) 

Jo Jr 

Vi(t) + V2 (t), (2.34) 

[ F m (t — s)dA(s) = — [ A(s)dF m (t — s), (2.35) 

(2.36) 

(2.37) 

(2.38) 


1 (s + Xj < t, \/j)dK(a(s),x), 

Z\{t) + z 2 (t), Zi(t ) := {Zij(t ),..., Z K ,i(t)), 7 = 1 , 2 , 

[ (Fk(t — s) — F m (t — s))dA(s) = [ A(s)d(F m (t - s) - F k (t - s)), 


o Jr 


(l(s + Xk < t) - l(s + Xj < t, Vj))dK(a(s),x) = -Mfc j2 (i) - F 2 (t)(2.39) 


The processes M 2 , Z 2 and V 2 are defined in the mean-square sense. This is in the same way 
as the limit process with respect to a standard Kiefer process for the G/GI/oo queue is defined in 
[28, 39]. The limit processes are characterized in the next subsection. 


2.4 Characterization of the Limit Processes 

In this section, we show the Gaussian property of the limiting processes (X,Y) and S when the 
arrival limit process is a Brownian motion. 

Theorem 2.4 (Gaussian Property). Under Assumptions 1 and 2, when the arrival limit process 
A is a Brownian motion, i.e., A{t) = c a B a (a(t )) for a standard Brownian motion B a , a positive 
constant c a > 0 and t > 0, the limiting processes (X,Y) and S in Theorem 2.3 are well-defined 
continuous Gaussian processes. For each t > 0, 

(X(t),Y(t)) = N(0,E(t)), and S(t) = N(0, a s (t)), 


where for j, k = 1,..., K, 


: = 

C OV (Xj(t), X k (t)) = 

J o s,t- s) + (c 2 a - 

1 )F?(t- 

-s)F fc c (t-s)] 

da(s), 

(2.40) 

:= 

Cov(Yj(t),Y k (t)) = 

1 " [(Fj,k(t ~ s,t- s) - F m (t 

-*)) 





+ ( c a 

~ !)( Fj(t - s) - F m {t - s))(F k (t - s) 

- F m (t - s)) 

] da(s), 

(2.41) 

4 y (0 :: 

= Cov{X j {t),%{t)) = 

= / (F k (t — s) — F jjk (t — s, 
Jo L 

t-s)) 





+ (c 2 a -l)(F^t-s)(F k (t~s] 

) - F m (t 

- s)))]do(s). 


(2.42) 

and 


fi 

rt 





a 5 ( t ) := Var(S(t )) 

= / F m (t - s)dd(s) + (cl - 
Jo 

1) J (. F m (t - s)) 2 da (, 

»)■ 

(2.43) 
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When the arrival rate function a(t) = A t for a positive constant A > 0, 

(X(t),Y(t)) => (X(oo),F(oo)) = IV(0,S(oo)) as t oo, 
lim f _1 Far(5(t)) = A c 2 , 

t—>oo 

where for j, k = 1, K, 

/•oo /*oo 

Ojfc(oo) : = A / s ) ds + A ( c a - !) / F j(s)Ff(s)ds, (2.44) 

JO JO 

4(oo) := A^°° [(F,- fc ( S , a ) - F m (a)) + (e£ - 1 )(Fj(s) - F m (a))(F fc (a) - F m (a))]da, (2.45) 

4 y (°°) : = A J o °° [(**M " + (4 - 1) (i?W(F fc (a) - F m (a))) ]da. (2.46) 

We make the following remarks on the Gaussian property of the limiting processes. 

(i) When we set = 1, the variance and covariance formulas coincide with those in the Poisson 
arrival case in Proposition 2.1. 

(ii) When K = 2 and c 2 a = 1, Cov(Yj(t),Yf.(t )) = 0, for t > 0 and k,j = with k / j , 

even if the service times of parallel tasks are correlated, since both terms inside the integral 
in (2.41) vanish. 

(iii) We emphasize the interesting structure of the variances of and Y^ and their covariances, 

k = 1, Recall that for G/GI/oo queues [28], the steady-state variance formula of the 

number of jobs in the system is given as the sum of two terms, the mean and the coefficient 
(c 2 — 1) multiplying an integral associated with the service time distribution; for example, 
when E[rji] < oo, the variance of the steady-state number of tasks in the k th service station 
is 

/•OO 

Var(X k (oo)) = \E[r 1 l] + X(c 2 a -l) ( Ff(s)) 2 ds , k = l,...,K. 

J o 

It turns out that the steady-state variance formula for the number of tasks in the waiting buffer 
for synchronization has the same structure; for instance, when E[rfA < oo for k = 1 
the variance of the steady-state waiting buffer size at the k th service station is 

/•OO 

Var(Y k (oo)) = X(E[7 1 l}-E[ V t}) + X(c 2 a -l) / (F^s) - Ff(s)) 2 ds, k = l,...,K. 

Jo 

The same structure also exists for the covariances between Xj and Y k , as shown in (2.42), for 
k. j = 1,..., K. 

(iv) The synchronized process does not have a Brownian motion limit, but its limiting process is 
Gaussian, and has the same variability as the arrival process when the arrival rate is constant. 
This can be also explained by regarding the synchronized process as the departure process 
of a G/GI/oo queue with the same arrival process and service times as the maximum of the 
service vectors (see [28, 39]). 
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To explore the impact of the correlation among the service times of each job’s parallel tasks 
on the system dynamics, we consider the case when the service vector rj l has the joint continuous 
distribution function 

F(x) = (1 - p) G(x k ) + pG (i min {xfcj'j (2.47) 

with a marginal continuous distribution function G(-), for 0 < p < 1, x k > 0 and k = 1 
Namely, the service times at the parallel stations have the same distribution, and are symmetrically 
correlated with a correlation parameter p £ [0,1). We state the mean and covariance functions of 
the performance measures studied above as functions of the parameter p in the following corollary. 

Corollary 2.1. Under the same assumptions in Theorem 2.4, when the service vector rf has the 
joint distribution function F in (2.47), for each t > 0 and k = 1,..., K, X k (t ) and V ar(X k (t)) are 
the same as in (2.23) and (2.40), respectively, 

Yk(t ) = (1 ~ p) [ [G(t - s)(l - (G(t - s))*- 1 )} da(s), 

Jo 

Var(Y k (t )) = J* [(1 - p)G(t - s)(l - (G(t - s ))^" 1 ) 

+ (1 - pf{c 2 a - l )(G{t - s)f (1 - (G(t - s)) 71 ” 1 ) 2 ] dd(s), 

Cov(X k (t),Y k (t )) = (c 2 a - 1)(1 - p) [ [G c (t - s)G(t - s)(l - (G(t - s))^ 1 )] dd(s), 

Jo 

for j, k = 1,..., K and j / k, 

Cov(Xj(t),X k (t)) = J* [(1 - p)(G c (t - s)) 2 + pG c (t -s) + (c 2 a - l)(G c (t - s)) 2 ]dd(s), 
C 0 v(Yj(t),Y k (t)) = J* [(1 - p)(G(t - s)) 2 ( 1 - (G(t - s)) K ~ 2 ) 

+ (1 - P) 2 (c 2 a - 1 )(G(t - s )) 2 (1 - (G(t - s)) K ^) 2 ]dd(s), 
G<w(Aj(i), y fc (t)) = (1 -p)J [G(t - s)G c (t - s ) 

+ (cl — 1 )G c (t — s)G(t — s) (l — (G(t — s))*- 1 ) j da(s), 

and 

S(t ) = J |^(1 - p)(G(t - s)) K + pG(t - s)jda(s), 

Var(S(t )) = f [(1 - p)(G(t - s)) K + pG(t - s)] da(s) 

Jo 

+ (c 2 - 1 ) / [(1 ~ p)(G(t-s)) R + pG(t,~ s)] 2 dd(s). 

Jo 
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We make several remarks on the impact of the correlation among the service vector. The mean 
and the variance of Xj~{t) are not affected by the correlation, but the covariances of Xj(t ) and X^(t) 
increase linearly in p for t > 0 and j, k = 1, ...K with j ^ k. The mean of Tfc(t) decreases linearly in 
p and the mean of S(t) increases linearly in p for t > 0 and k = 1, K. The covariances of Yj(t) and 
Yk(t) decrease in p. in the order of (1 — pj 2 , but the covariances of Xj(t) and !&( t) decrease linearly 
in p for t > 0 and j,k = 1, ...K . The variance of S(t) increases in p, in the order of p 2 , for t > 0. 
The intuitive interpretation for these observations is that positive correlation makes the parallel 
tasks more likely to finish close to each other so that the waiting time for synchronization becomes 
less and more jobs are synchronized. It is also important to emphasize that the covariances of Yj[t) 
and Yk(t) and the covariances of Xj(t) and Y^(t) decrease in different orders in the correlation 
parameter p for t > 0 and j,k = 1, ...K. The same observations hold for the associated steady-state 
performance measures. 


2.5 Comparison with a fork-join network with ES 


We make a comparison with an associated fork-join network with ES. We use superscript “ES” in 
the corresponding processes for this model. Let the arrival and service processes be the same as 
the model described above. The only difference is the synchronization constraint. Here tasks are 
not tagged with a particular job, so that whenever there are tasks completed at all parallel service 
stations, the oldest completed task at each waiting buffer for synchronization will be synchronized. 
It is evident that when the arrival process A(t) is Poisson, the processes Y k {t) and S ES (t) do not 
have a Poisson distribution at each time t > 0, k = 1 In this case, for each k = 1 

X' k ,ES and ]J™' ES w iii have the same representations as in (2.6) and (2.9), but the processes S n,ES 
and y£’ ES become 

S n ’ ES (t) = min W;’ ES (t)}, t> 0, (2.48) 

1<J<K J 

and 

Y?' ES (t) = D n k ' ES (t) - S n ' ES (t) = D n k ’ ES {t) - min {D]’ ES (t)} : t > 0. (2.49) 

Thus, at any time, one of the waiting buffers for synchronization should be empty. It is evident that 
the processes S n ' ES and Y k ' E cannot be represented as a single integral of the multiparameter 
sequential empirical process K n as in equations (2.15) and (2.14), respectively. 

We now discuss more on the comparison for the steady-state mean values of the fluid limits of 
these processes when the arrival rate is constant. In the ES model, the synchronization process 
S ES can be represented as the minimum of the departure processes from all parallel stations, and 
these departure processes are dependent due to the correlation of service vector of each job. Thus, 
we are unable to obtain a distributional approximation of the processes S ES and Y ES , k = 1,..., K. 
However, for each t > 0, by applying the previous results on G/GI/oo queues [28], we can obtain 
the mean values of the fluid limit Y ES (t), k = 1,..., K, and S ES (t ): 


Y ES (t) := A 


Fu(s)ds — min 
1 <j<K 


Fj(s)ds 


Y, 


f S (oo] 


■■= A ( 


. max^lE^]} - E[ri k ] 


as t —> oo, 


(2.50) 
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S ES (t) := A min / f FTsW.sl = Xt — X max \ f F c (.s)dsl , and lim — — 1 — = A. 

(2.51) 

Recall that the steady-state mean value of the waiting buffer for synchronization in our model 
F fc (oo) = A (E[nl] - £7[T7^]) in (2.27), k = 1 denoted as Y EES ( oo) for the comparison 

purpose. It is evident that the average waiting buffer sizes for synchronization under NES constraint 
are larger than those under ES constraint, even though the total synchronization throughput rates 
are the same, lim^oo S ES (t)/t = lim^oo S NES (t)/t = A. We also observe that when the parallel 
service times are perfectly positively correlated, the difference Y^ ES ( oo) — Y ES ( oo) becomes zero 
for k = 1, We summarize this comparison result in the following proposition. 

Proposition 2.2. Under Assumptions 1 and 2, when a(t) = Xt for a positive arrival rate A > 0 
and E[r]jf\ < oo for k = 1,..., K , 

Yk ES {°°) ~ Y ES (oo) = X(E[r,l] - max K {E[ V ]]}) > 0, for k = 1,..., K. (2.52) 

By the extreme value theory, if the service vector has i.i.d. components such that the service 
time distribution lies in the domain of attraction for Gumbel extremal distribution, then we have 
a K(.Vm ~ frftr) =>• Z as K —> oo, where Z has a Gumbel distribution, and ax and bx are constants 
depending on K ; see Chapter 1 in [31]. The Gumbel distribution has cdf P(Z < z) = e~ e “, z > 0, 
with mean E[Z] = 7 ~ 0.5772, the Euler-Mascheroni constant, and variance Var(Z) = 7 r /\/6 ~ 
1.2825. For one example, if the service vector has i.i.d. components of an exponential distribution 
with rate 1, then ax = 1 and bx = In K (see Example 1.7.2 of [31]), for k = 1,..., K, 

Y k NES ( 00 ) - Y es ( 00 ) = A ^ l - 1^ « A(ln (K) - 1) as K ^ 00 . (2.53) 

For another example, if the service vector has i.i.d. components of a lognormal distribution 
LN(0 ,1), we have, for k = 1,..., K, 

Y k rES (oo) — Y es (oo) « X^/ax + bx — e 1 / 2 ) as K —>• 00 , (2-54) 

where a/v and are (see Example 1.7.4 of [31]): 

ax = (2 In 77 ) 11/2 exp | — (2 In A") 1 / 2 + 0.5(21n/t ) _1 / 2 (lnlnRT + ln(47r))| , 

and 

bx = exp |(21nil') 1 / 2 — 0.5(21n77) _1 / 2 (lnlniA" + ln(47r))| . 

2.6 Numerical Example 

In this section, we provide a numerical example with two parallel tasks (K = 2), comparing our 
approximations with simulations. We let the arrival process be renewal with arrival rate A = 100 
and the SCV c 2 = 5. The service times of the two parallel tasks are assumed to be a bivariate 
Marshall-Olkin hyperexponential distribution, which is a mixture of two independent bivariate 
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Marshall-Olkin exponential distributions [34]. A bivariate Marshall-Olkin exponential distribution 
function Fmo{x,u) for the random vector (X, Y) can be written as F^ IO {x 1 y) := P(X > x,Y > 
y) = exp(— /i\x — p 2 y — pi 2 (x V y)), x, y > 0, where three parameters p 2l p\ 2 are such that the 
two marginals are exponential with rates p± + p\ 2 and p 2 + y\ 2 and their correlation p = pi 2 /(pi + 
P 2 + P 12 ) £ [0,1]- We denote MO(X\, X 2l p) for a bivariate Marshall-Olkin exponential distribution, 
where Ai and X 2 are the rates for the marginals, and p is the correlation parameter, for which the 
parameters p\ = (Ai - p\ 2 )/{l+p), p 2 = (\ 2 - p\±) / (1 + p) and p\ 2 = (p(Ai + A 2 ))/(l + p). In the 
numerical example, we take a mixture of MO(4/5,1, p\) with probability 0.4 and M 0(6/5,6/5, p 2 ) 
with probability 0.6, such that the means of the two hyperexponential marginals are m s , 1 = 1 and 
m s ,2 = 0.9. By setting p\ = p 2 = 0, we have two independent parallel service times, and by setting 
pi = 0.7 and p 2 = 172/679, we obtain that the correlation (see the correlation formula in §5.2 [40]) 
between the two parallel service times is equal to 0.5. 

In Table 1, we show the approximation values for the mean, variance and covariance of Xk 
and I*,, for k = 1,2, and compare them with the corresponding simulated values. To estimate 
the simulated values, we simulated the system up to time 40 with 4000 independent replications 
starting with an empty system, which we call one experiment. In each replication, we collected data 
over the time interval [20,40] and formed the time average (the system tends to be in steady state 
in less than 5 time units). We conducted 5 independent experiments and took sample averages as 
estimations for simulated values. To construct the 95% confidence interval (Cl), we used Student 
t-distribution with four degrees of freedom. The halfwidth of the 95% Cl is 2.776s5/\/5, where S 5 
is the sample deviation. 

We make several remarks for the numerical example. First, our approximations match very well 
with the simulated values. Second, the size of waiting buffers for synchronization is quite large, of 
the same order as the number of tasks in the service stations. Third, we find that when the two 
parallel tasks are positively correlated, the mean and the variance of X^s are the same as those 
in the independent case, while the covariance between X\ and X 2 gets larger, the mean and the 
variance and covariances of Yj/s and the covariances between X^ and Yj become smaller than those 
in the independent case, j,k = 1,2. These are also consistent with the observations in Corollary 
2.1. Note that this numerical example is more general than that considered in Corollary 2.1. 

3 The multi-server fork-join network model 

3.1 Model and Assumptions 

In this section, we present a detailed description of our multi-server fork-join network model. We 
consider a fork-join network with a single class of jobs, and each job is forked into K ( I\ > 1) 
parallel tasks. Each task is processed in a service station with finite servers under the non-idling 
FCFS discipline. Namely, a newly arriving task immediately gets served if there is an idle server in 
that station, and joins the back of the queue otherwise, and the task waiting for the longest in the 
queue enters service as soon as a server in that station becomes available. After service completion, 
each task will join a waiting buffer for synchronization associated with each service station, and 
when all tasks of the same job are completed, they will be synchronized and leave the system. Here 
we assume that the synchronization process takes zero amount of time. 

Let A := {A(t) : t > 0} be the arrival process of jobs after time 0. Let t* be the arrival time of 
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Table 1: Comparing approximations with simulations in a stationary model 


(Xi,x 2 ) 

(£[Xi],£[X 2 ]) 

(Var(Xi), Var(X 2 )) 

Cov(X i,X 2 ) 

P = 0 

Sim. (95% CL) 

(99.99 ± 0.17 , 89.98 ± 0.12) 

(296.26 ± 0.66, 269.46 ± 0.70) 

234.14 ± 0.66 

Approx. 

(100.00, 90.00) 

(296.00, 269.27) 

233.99 

p = 0.5 

Sim. (95% CL) 

(99.98 ± 0.04, 89.99 ± 0.04) 

(296.08 ± 0.57, 269.23 ± 0.80) 

256.34 ± 0.43 

Approx. 

(100.00, 90.00) 

(296.00, 269.27) 

256.30 


(Yi ,y 2 ) 

(ElYrfElY,}) 

(Var(Yi), Var{Y 2 )) 

cw(yi,y 2 ) 

P = o 

Sim. (95% CI.) 

(43.18 ± 0.05 , 53.20 ± 0.10) 

(70.12 ± 0.20, 89.85 ± 0.40) 

31.53 ± 0.30 

Approx. 

(43.20, 53.20) 

(70.31, 90.08) 

31.55 

p= 0.5 

Sim. (95% CI.) 

(20.89 ± 0.01, 30.88 ± 0.02) 

(27.14 ± 0.15, 42.23 ± 0.35) 

8.36 ± 0.07 

Approx. 

(20.89, 30.89) 

(27.05, 42.23) 

8.31 


(X,Y) 

Ccw(Ai, Yi) 

Cov{ A'i,y 2 ) 

Cov(X 2 ,Y 1 ) 

Cov(X 2 ,Y 2 ) 

P = o 

Sim. (95% CI.) 

60.80 (± 0.59) 

122.87 (± 0.61) 

99.21 (± 0.42) 

64.56 (± 0.54) 

Approx. 

61.09 

123.10 

99.85 

64.57 

p = 0.5 

Sim. (95% CI.) 

28.72 (± 0.33) 

68.37 (± 0.73) 

47.51 (± 0.42) 

34.49 (± 0.44) 

Approx. 

28.67 

68.37 

47.41 

34.44 


the i th job, i € N, that is, Aft) = max{i > 1 : 77 < t} for t > 0 and A(0) = 0. Let Nk be the number 
of servers at service station k, k = 1, I\. Each job brings in a /L-dimensional service vector, 
representing the service time at each service station, which can be correlated. Let rf := (rj\ ,... ,rf K ) 
be the service vector of the job that arrives at time r*, i 6 N, where r/t is the service time at the 
k th service station. We assume that the sequence {rf : i > 1} is i.i.d., and let the joint distribution 
function of rf be F(x) = F{x\,..., xk ) for x^ > 0, k = 1, Let F c (x) := P(rfi > x\ , ...,rf K > 

xk), for x \,...,xk > 0. Their marginal distributions are Tfc(-) with mean 1/Hk £ (0, oo), for 
k = 1,..., K. Let rf m := maxjbyj,..., rf K } and F m (x ) := P{rf m < x) = P(rfj < x,Vj) = F(x,...,x ) 
for x > 0. (Throughout this paper, we use “m” to index quantities and processes associated with 
the maximum.) We make a regularity assumption on the service time distributions for the parallel 
tasks. 

Assumption 1. The joint distribution function F(x) of the service time vector rf, i G N, is 
continuous. 

State Descriptors. Let := (Xfc(f) : t. > 0} be the process counting the number of tasks at 
the service station k, and Yj, := {Yfe(t) : t. > 0} be the process counting the number of tasks in 
the waiting buffer for synchronization (unsynchronized queue) after service completion at service 
station k, k = 1 Denote X := (X \,..., Xk) and Y := (Y\, ...,Yk). Let S := { S(t ) : t > 0} 

be the process counting the number of synchronized jobs by each time t > 0. In addition, let 
Qk '■= { Qk{t ) : t > 0} and B^ := {B^ft) : k > 0} be the processes representing the queue length 
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and the number of tasks in service at station k, respectively, k = 1, K. Let D := {D^ft) : t > 0} 
be the cumulative service completion (departure) process at service station k, k = 1, K. Denote 
Q ■■= (Qi, Qk), B := (Sj, ...,B k ), and D := (D 1; D K ). 

A Sequence of Systems. We consider a sequence of the above fork-join networks, indexed by 
superscript n and let n oo. We assume that each service station is operating in the many-server 
heavy-traffic asymptotic regimes, where the arrival rate of jobs and the number of servers get large 
appropriately while the service time distributions are fixed. In establishing the FLLN, we allow 
the arrival rate to be time-dependent. In establishing the FCLT, we will assume that each service 
station is operating in the Haffin-Whitt (QED) regime, so that it is critically loaded with a constant 
arrival rate (see Assumption 4 for the precise definition). For any process A, we use X n to represent 


the associated process in the sequence of the fork-join networks. 

Some Fundamental Flow Balance Equations. For each service station k. k = 1 and for 

each t > 0, we have the following flow conservation equations: 

x%(t) = m)+Qm, (3.i) 

XZ(t)=X%(0) + A n (t)-mt)> (3-2) 

Y?{t) = Y?(0) + DW)-S n {t). (3.3) 

The non-idling condition implies that for each k = 1,..., K and t > 0, 

Bi(t) = X$(t)ANi, Ql(t) = (XZ(t)-NZ)+. (3.4) 


In addition, we have the following flow balance equation, for each k, k! = 1,..., K , k A k ', and t > 0, 


jq?(t) + i?(t) = *S(t) + im (3-5) 

that is, the total numbers of tasks in each service station and its associated waiting buffer for 
synchronization are equal at all time, and are equal to the total number of jobs in the system. 


3.2 Fluid Limit 

In this section, we present the fluid limit for the fork-join network. We assume that the system 
starts from empty and allow the arrival rate to be time-dependent. 

Assumption 2. There exists a continuous nondecreasing deterministic real-valued function d(t) 
on [0, oo) with a(0) = 0 such that 

A n (t ) := n~ 1 A n (t) =>- aft ) in D as n —> oo. (3-6) 


We also make the following assumption on the numbers of servers. 

Assumption 3. For k = 1,..., K, NJf := NJf/n —> Nk > 0 as n —> oo. 

Under the empty initial condition, we can write the processes X'f(t), YA(t), k 
S n (t ) as 


A n (t) 


x k(t) = 1 ( r i l + w k’ 1 + 4 > t), t> o, 


i=1 


1, ...,K, and 

(3.7) 
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(3.8) 


A n (t) 

Y k (t) = ^ 1 (rf + + ??[. < t, rf + + jyjb > f, for some k' / fc), t> 0, 

1=1 

A n (t) 

S n (t) = ^ l(r” + w;]?’* + 7/j{. < t, V&: = 1,..., IF), f > 0, (3.9) 

i=l 

where vj?’ 1 is the waiting time of the i th arrival at station k , i £ N. 

In addition, for k = 1 let E%(t) be the number of tasks that have entered service at 

station k by time t, t > 0, and set := {E^(t) : t > 0}. Denote E n := (E 1 ”,..., E^). For each 
service station k = 1, K, we also have the balance equation 

m) = A n (t) - QW) = A n (t) - {xm - NJ!) + , t > 0. (3.10) 

Define the fluid-scaled processes X n := n~ l X n for X n = X n ,Y n ,S n ,E n ,Q n ,B n ,D n . We now 
state the FLLN for the fluid-scaled processes. 

Theorem 3.1. Under Assumptions 1-3, 

(. A n ,X n ,Y n ,S n ,E n ,Q n ,B n ,D n ) => ( a,X,Y,S,E,Q.B,D ) (3.11) 

in p6A'+2 

as n oo, where the limits are all deterministic functions: a is the limit in (3.6), 
(. E,X,Y,S ) is the unique solution to the following: fort > 0 and k = 1 

X k (t) = T Ff(t - s)da(s) + f(X k (t - a) - N k )+dF k (s), (3.12) 

do Jo 

F k (t) = d(t) - (X k (t) - N k )+, (3.13) 

S(t) = [...[( min {E k (t - s k )}] dF(si,...,s K ), (3-14) 

J0 Jo \k=l,-,K J 

Y k (t) = f F k (t - s)da(s) - [\x k (t - s) - N k ) + dF k (s) - S(t), (3.15) 

Jo Jo 

and the limits Q, B and D satisfy 

Q k (t) = (X k (t) - N k ) + , B k (t) = X k (t) A N k , D k (t) = a(t) - X k (t). (3.16) 

It is easy to check that for each k = 1,..., K, the limit X k (t) also satisfies the following equation: 

X k (t) = a(t ) - [ E k (t — s)dF k (s), t> 0. (3-17) 

Jo 

When a(t) = [q \(s)ds and the service times are exponential (independent or dependent), where 
A(-) is a positive function, for each k = 1, the fluid limit X k in (3.12) and (3.17) becomes an 
ordinary differential equation (ODE) [38], but the fluid limit Y k in (3.15) does not have an ODE 
representation. We remark that the fluid limit X k for each k = 1,..., K depends only on the marginal 
distribution F k , while the fluid limits Y k , k = 1,..., K, and S depend on the joint distribution F. 
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However, as the FCLT (Theorem 3.2) below shows, the limits for all these processes in the diffusion 
scale will depend on the joint distribution F. 

When the arrival rate is constant and each service station is underloaded or critically loaded, we 
give a corollary on the steady states of the fluid limits. The proof follows from a direct calculation 
and is omitted. It is evident that correlation among service times of parallel tasks only affects the 
steady state of Y but not that of X. 

Corollary 3.1. Under Assumptions 1-3, if the arrival rate is constant, d(t) = A t, for A satisfying 
0 < A < N k p k for all k = 1,..., K, 

{X(t),Y(t),Q(t),B(t)) —> (X(oo),Y(oo),Q(oo),B(oo)) as t — > oo, 

and 

~(D(t), E(t), S(t)) —>• A := (A,..., A) as t — > oo, 

where 

X k (oo) = B k { oo) = A E[ril] = X /Uk, %{ oo) = A [E[p} n ] - E[r)\\), Q k {oo) = 0. 


3.2.1 Numerical Examples 

We give two numerical examples to show the effectiveness of fluid approximations comparing with 
simulations, when I\ = 2. We let the arrival process be Poisson with time-varying rate A (t) = 
200 + 120 sin(t), t > 0. The numbers of servers in stations 1 and 2 are N] = 300 and IV 2 = 340, 
respectively. In the first numerical example, the service times of the two parallel tasks are assumed 
to have a bivariate Marshall-Olkin exponential distribution [34], For our first numerical example, 
we set the service times to be MO(l,0.9,p) such that the service times of the two parallel tasks 
have exponential marginals with means 1 and 10/9 in stations 1 and 2, respectively, and their 
correlation is p. The numerical results with p = 0 and p = 0.5 are provided in Figure 3(a), marked 
with “ind.” and “corr.”, respectively. In the second numerical example, we let the service times of 
the two parallel tasks have a bivariate Marshall-Olkin hyperexponential distribution [40], which is 
a mixture of two independent bivariate Marshall-Olkin exponential distributions. Specifically, we 
take a mixture of MO{ 4/5, l,p\) with probability 0.4 and MO(6/5, 27/32, P 2 ) with probability 0.6, 
such that the two parallel service times have hyperexponential marginals with the same means as 
the first example. By setting pi = p -2 = 0, we have two independent parallel service times, and 
by setting p\ = 0.7 and P 2 = 521/1232, we get the correlation between the two parallel service 
times to be 0.5. In Figure 3(b), we show the numerical results with p = 0 (“ind.”) and p = 0.5 
(“corr.”). To calculate the simulated values, we simulated the system up to time 20 with 500 
independent replications starting with an empty system. We make two remarks from numerical 
results. First, the fluid approximations match very well with the simulated results. Second, the 
positive correlation among parallel service times does not affect X k , but reduces Y k , for k = 1,2. 

3.3 FCLT in the Halfin-Whitt regime 

In this section, we study the fork-join network with NES in the Halfin-Whitt regime, which requires 
that each service station operates in a critically loaded regime asymptotically. Specifically, we 
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Figure 3: Comparison of fluid approximations with simulations. 

assume the following. Let \ n be the arrival rate of jobs such that A n := A n /n —> A > 0 as n —> oo, 
and set NJ} := nN k , where N k e N, and p k := \ n /(p k N k ) for each k = 1, ..., A'. 

Assumption 4. For each k = 1,..., A', A = N k pk and, y/n(l — p k ) —>• > 0, as n —>• oo. 

The arrival processes A n = {A n (f) : t > 0} satisfy an FCLT. 

Assumption 5. There exists a stochastic process A with continuous sample paths satisfying 

A n (t) — \ n t 

A n (t) :=- -j= -=4- Aft) in D as n —> oo. (3.18) 

\Jn 

It follows from (3.18) that we have the associated FLLN: 

A n (t) =>• Xt in D as n —> oo. (3.19) 

We now describe the non-empty initial conditions. Due to the complexity from initial conditions, 
we focus on the case of K = 2, but our approach can be extended to K > 2. For convenience, we 
use the notation k' to denote its counterpart, i.e., k' = 1 (k' = 2, respectively) if k = 2 (k = 1, 
respectively), for k = 1, 2. At time 0—, there are A^(0) tasks at service station fc, and Y£( 0) tasks 
in its associated waiting buffer for synchronization, for k = 1,2. Let A n (0) := (X™(0), Xhf (0)) and 
F n (0) := (T 1 n (0),T 2 n (0))- Recall the flow balance equation (3.5). At time 0—, 

** (0) + Y k n ( 0) = X%(0) + Y£{ 0), k = 1, 2, (3.20) 

which is equal to the number of jobs in the system. Note that AT(0) > 171(0) for each k = 1, 2, 
since tasks in the waiting buffer associated with station k' for synchronization must be in station 
k, either in service or in queue. Let B k { 0) := min(X^(0), Njf) and <5^(0) := (Xj?(0) — NJf) + be the 
number of tasks in service (busy servers) and the queue length at station k at time 0—, respectively, 
k = 1,2. We also assume that 171(0) < B k ( 0) for k = 1,2. This is not a restrictive assumption, 
because in the Halfin-Whitt regime, waiting times for service at each station are 0(1/y/n) and 
service times are 0(1), and jobs that have completed tasks in one station and joined its waiting 
buffer for synchronization have their associated tasks receiving service in the other station with 
probability one asymptotically. 
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Let J n ( 0) := min/ i . = i ! 2{-B^(0) — 17j(0)} be the number of jobs whose both tasks are in service 
at time 0—. Then Z k ( 0) := B k ( 0) — Y$( 0) — J n (0) represents the number of jobs in the system at 
time 0— whose task k is in service but whose task k' is in queue waiting for service, k = 1,2. Let 
I n ( 0) := Q?(0)AQ2(0) be the number of jobs (if any) whose both tasks are in queue at their service 
stations at time 0—. Then R k (0) := <5)1(0) — I n ( 0) represents the number of jobs (if any) whose 
task k is waiting in queue for service while whose task k' is in service, k = 1,2. (Note that our 
assumption above implies that if a job is waiting in queue at station k, its parallel task can be either 
in queue or in service at station k'.) By our definition, we can see that Z k { 0) = R k ,( 0), fc = 1,2. Set 
R n ( 0) := (i?”(0), i ?2 (0)) and Z n ( 0) := (Z”(0), Z^ (0)). We also obtain a decomposition for X%(0): 

xm = BU 0) + QUO) = Y£( 0) + J n (0) + ZU 0) + I n ( 0) + R n k { 0), k = 1, 2. (3.21) 

We let {u>U : i = 1, Q k (0)} be the sequence of remaining waiting times of the tasks in 

station k at time 0—, k = 1,2. It is in the order of their positions in queue: uu ’ 1 is the remaining 
waiting time of the task in the front of the queue while is that for the task in the end 

of the queue at station k at time 0—, k = 1,2. Let {fj l k : i = 1,££(0)} be the sequence 

of remaining service times of the tasks in station k at time 0—, for k = 1,2. Let {rf k : i = 

1 , ...,< 5 )!( 0 )} be the sequence of service times of the tasks in station k that are in queue at time 
0—, k = 1,2. Without abuse of notation, we use {fj k Yk : i = 1,..., YJ!?(0)}, {rj k J : i = 1,..., J n (0)} 
and {fj k Z : i = 1,..., Z k ( 0)}, which are partitioning subsets of { ff k : i = 1,..., Bjjj(0)}, to represent 
the remaining service times of the tasks in station k at time 0 — corresponding to the quantities 
Y k ,{ 0), J n (0) and Z k ( 0), respectively, k = 1,2. Similarly, we use {wU’ 1 ■ * = l,---,-f n (0)} and 
{w k ,l,R : i = 1 , •••, -Rfc(O)}, which are partitioning subsets of {w k l : i = 1 ,..., Qjj( 0 )}, to represent the 
remaining waiting times of the tasks in station k at time 0 — corresponding to the quantities I n ( 0 ) 
and R k ( 0), respectively, k = 1,2. Finally, we use {if k : i = 1,..., / n (0)} and {rf k R : i = 1,..., i?^(0)}, 
which are partitioning subsets of {rf k : i = 1 ,..., < 3 ^( 0 )}, to represent the service times of the tasks 
in station k corresponding to the quantities I n ( 0) and R k ( 0) in queue at time 0—, respectively, 
k = 1,2. We assume that these initial quantities are independent of the arrival process A n and the 
service times of new arrivals after time 0 . 

We can now give a representation for the processes X n , Y n and S n : for t > 0 and k = 1,2, 


xm 


S n (t) 


and 


v 2 n ( 0) _ n n ( 0) J n ( o) 

<t)+ Y 1 ^ 2 - *) + Y 1 ^i’ J - 

i =1 i =1 i =1 

z?( 0 ) zg( 0 ) 

+ < t, W?/ Ji + <t)+ Y mY' R + nf < t, kif 

i =1 i =1 

I n i 0) A"(t) 

+ y my 4 +Vj 1 <tyj)+ Y + w 7 l + v) ^ 

i —1 2—1 


(3.23) 

<t) 


Y k m = Y k uo)+xm)+Am-xm-sm- ( 3 . 24 ) 
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We use the convention that X^=i = 0 throughout the paper. 

We impose the following assumptions on the initial quantities. 

Assumption 6. There exists (Yi(0), Y^O)) £ such that 

(A n (0),y"(0)) :=n _1 (JC n (0),y n (0)) =*► (X(0),F(0)) in M 4 as n -> oo, 

whereX( 0) := (Ni,N 2 ) andY( 0) := (Yi(0), 7^(0)). There exist random vectors X (0) := (Xi(0), ^ 2 ( 0 )) G 
M 2 andY( 0) := (Yi(0), ^ 2 ( 0 )) G M 2 such that 

(X n (0),y n (0)) := Vn{X n (0) -X(0),y n (0) -y(0)) =► (X(0),y(0)) in M 4 as n 00. 

This assumption implies that the associated fluid-scaled initial quantities 

(J n (0), (0), / n (0), R n (0)) := n^ 1 (J n ( 0 ), Z n ( 0 ), I n ( 0 ), R n ( 0 )) => (J( 0 ), Z(0), 7(0), 72(0)) 

in M 6 as n —> 00, where 


J(0) := Ai-y 2 (0) = A 2 -yi(0), 2 ( 0 ) := (Z 1 (0),Z 2 (0)) := (0,0), 7 ( 0 ) := 0, R( 0) := (0,0). 

Define the associated diffusion-scaled quantities (J n (0), 2 " (0), 7 n (0), R'\o)) by 


> ( o ):= J "(°>-" J < 0 > , ^ (0):= ^m 


r( 0):=ih£l, flj(0) := k = 1,2. 


Then Assumption 6 implies that 

^J"(0),Z n (0),7 n (0),^"(0)) (j(0),Z(0),7(0),.R(0)) in M 6 as n ^ oo, 

where 

7(0) := min{-(X fc (0))- - ^(0)}, Z k ( 0) := -(X fc (0))“ - Y k ,(0) - 7(0), k = 1,2, 

k= 1,2 

7(0) := min (A fc (0))+, R k ( 0) := (X fc (0))+ - 7(0), 7 = 1,2. 

k= 1,2 

Let 

1 7* 

7i-,eW : = 1=TTT / F k( s ) ds ’ * ^ °- 
7o 

be the equilibrium distribution associated with F k , k = 1,2. 

Assumption 7. Fork = 1,2, {;/)(. : i G N} is a sequence ofi.i.d. random variables with distribution 
F ke and for each i G N, j)) and rj | are independent. {rf k : i G N} is a sequence of i.i.d. random 
variables with distribution F\ for each i G N and k = 1,2. {(rff 1 , rff 1 ) : i G N} is a sequence of i.i.d. 
random vectors with a joint distribution F(-,-). {(rff ,rj k , ) : i G N} is a sequence of i.i.d. random 
vectors with independent components, k = 1,2. 

Finally, we also make an assumption for the residual waiting times {vo^' 1 : i = 1,..., Q k (0)}, 
k = 1, 2. 
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Assumption 8. The residual waiting times of the tasks in queue {wf’ 1 : i = 1, Qf( 0)} ; k = 1, 2, 
converge to zero a.s. as n —> oo. 

We define the diffusion-scaled processes X" := (X", Xf), Y" := (F”, Kf) and S n by 


xm := 


xm~ N k 


Y k n {t) : = 


W)-W) 


5 n (f) : = 


5 n (f) - 5"(t) 


for k = 1,2, where 

S n (t) := nS°(t)+ \ n f f ((t - si) A (t - s 2 )) dF(s\, s 2 ), (3.26) 

Jo Jo 

S°(t) := Y 2 (0)F lie (t) + Fi(0 )F 2 , e (t) + J(0)F 1>e (f)F 2>e (t), (3.27) 

Y£(t) := nYfc(O) + A n t - S n (t). (3.28) 

From the balance equation for Yff in (3.24), we can rewrite Yff as 

Y£(t) = YfJ (0) + X£(0) + A n (t) - Xjf(t) - S n (t), t > 0, k = 1, 2. (3.29) 

Recall EV(t) is defined as the cumulative number of tasks entering service by time t > 0 at 
station k, k = 1, 2, assuming the system starts empty in §3.2. Without abuse of notation, in §3.3 
related to the FCLT, we let E%(t) be the number of new arrivals after time 0 whose task k has 
entered service by time t > 0 at station k, k = 1,2. 

Define the diffusion-scaled processes (E ,Q ,B ,D ), E := (Ef, Elf ), Q := (Qf, Qf), B := 
and I)' 1 := (D?,D%), by 

z^n/j-A _ \ri-L 

m) ■■= - jA2 7 =—, Qm ■■= (xm)) + , m) ,= -(. xmr , 

\Jn 


L>m := Xjf(0) + A n (t) - X%(t), t> 0, k = 1, 2. (3.30) 

For si, s 2 > 0, let 

£ n (si,s 2 ) := -i((fi; l ( Sl )A^ 2 ))-A’*( Sl A S2 )) 

V n 

= (F'i'(si) + (A"/\/n)( s i — si A s 2 )) A (Elf (s 2 ) + (A n / \fn)(s 2 — s\ f\ s 2 )).(3.31) 

Before we present the FCLT for the fork-join network with NES in the Halfin-Whitt regime, 
we provide some preliminaries for the limit processes. The limit processes will be functionals of 
a generalized multiparameter Kiefer process, as a limit of the multiparameter sequential empirical 
process driven by the service time vectors of new arrivals. Define the multiparameter sequential 
empirical processes K n := {K n (ti, t 2 , x) : t\ > 0, f 2 > 0, x G M^_} by 

\nti\ A \r1t2} 

K n (t 1 ,t 2 ,x) :=-= V ( 1 ( 17 * < x) - F(x)). (3.32) 

Vn 

We prove the convergence of K n in the space D([0, oo) 2 ,D([0, oo) 2 , M)) endowed with a generalized 
Skorohod Ji topology defined in [18] in Proposition 3.1. 
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Proposition 3.1. Under Assumption 1, 

K n (ti,t 2 ,x) =>• K(ti,t 2 ,x) in D([0, oo) 2 , D([0, oo) 2 , M)) as n —> oo, (3.33) 

where K(ti,t 2 ,x) is a continuous Gaussian random field, called a generalized multiparameter Kiefer 
process, with mean E[K(t\, t 2 ,x)\ = 0 and covariance function 

Cov(K(si, S 2 ,x),k(ti,t 2 ,y)) = (si A s 2 A t\ A t 2 )(F(x Ay) - F(x)F(y)), (3.34) 

for Sfc, tk > 0, k = 1, 2, and x,y G M+. 

We define the processes W k := (Wfc(t) : i > 0}, := {Wf(t) : t > 0} and W := {W(t) : t > 0} 

as integral functionals of K: for t > 0, k = 1, 2, 

Wfc(t) := f f f l(sfc + .Tfc < t)dk(Xsi, As 2 ,x), (3.35) 

Jo Jo J«.\ 


W(t) := / / / l(sj + Xj <t,'ij)dk(Xsi,Xs 2 ,x), 


(3.36) 


and 


Wf(f)-.= W k {t)-W(t)= / / l(s fc + x k < t, s k > +x k r > t)dk{\si,Xs 2 ,x), (3.37) 


where the integrals are defined in the sense of mean-square limits (see the precise definition in §??). 

Proposition 3.2. The processes W k , Wf and W are well-defined continuous Gaussian processes 
with mean zero, and for 0 < s < t and k = 1,2, 

E[(W k (f) - W fc (s)) 2 ] = A F k (t -u)- F k (s - «))(1 - F k (t - u) + F k (s - u))du , 

Jo 


E[(W(t) - W(s)) 2 ] 


A [ [ [AF((s - si,s - s 2 ); (t - si,t - s 2 ))] 
Jo Jo 


x [1 - A F((s — si,s — S 2 ); (t — si,t — S 2 ))]d(-Si A s 2 ), 


(3.38) 


E[(Wf(t) - Wf(s)) 2 ] = E[(W k (t) - W k (s)f ] + E[(W(t) - W(s)) 2 ] 

- 2A / / [F(t - si,t - s 2 ) - F k ,k’(s - s k ,t - s k >) 

Jo Jo 

+ ( F k (t - s k ) - F k (s - s k ))(F(s - si, s - s 2 ) - F(t - si,t - s 2 ))]d(si A s 2 ), 
and covariance functions 


Cov(W k (t),W k '(t )) 


A [ [ [F(t - si,t - s 2 ) - F k (t - s k )F k fit - s k >)]d(si A s 2 ), 
Jo Jo 
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Cov(W k {t),W k ,(t)) = A [ f [F k (t - s k )F(t -si,t- s 2 ) - F k (t - s k )F k >(t - s k >)]d(si A s 2 ), 

Jo Jo 

Cov(W k (t), W(t)) = A [ [ [F(t - si,t - s 2 ) - F k (t - s k )F(t - si,t - s 2 )]d(-si A s 2 ), 

Jo Jo 

Cov(W k (t),W(t)) = A f [ [(F(t - si,t - s 2 )f - F k (t - s k )F(t - si,f - s 2 )]d(si A s 2 ), 

Jo Jo 

where F k>k t(x,y ) := P{rf k < x, rf k , < y) for x, y G M+, and 

AF(x;y) := F(y 1 ,y 2 ) - F(xi,y 2 ) - F{y u x 2 ) + F(xi,x 2 ), x,y <E K 2 , a; < y. 

In addition, let 1/ := {U(t) : t £ M+} be a continuous two-parameter Gaussian process with 
mean zero and covariance function: 

Cov(U(s),U(t )) = (Fi, e (si A ti)F 2>e (s 2 A t 2 ) - F lje {s{)F 2te (s 2 )F^ e {t{)F 2te (t 2 )), (3.39) 

for s := (s\,s 2 ) £ and t := (fi,f 2 ) G M+. Define U k := (74(1) : t > 0}, for k = 1,2, by 

74(f) := Z7(f,oo), U 2 (t) := U(oo,t), t> 0, (3.40) 

and without abuse of notation, we denote 17(f) = U(t,t), t > 0. Note that the processes 114, ID): 
and ID are independent with 17, as well as 74, k = 1, 2. 

We are now ready to state the FCLT. 

Theorem 3.2. Under Assumptions 1 and 4-8, 

(. A n ,X n ,Y n , S n ,E n ,Q n ,B n ,b n ) => (A,X,Y,S,E,Q,B,d') (3.41) 

in D 14 as n —>• oo, where A is in (3.18), X, Y and S are the unique solutions to the following set 
of stochastic integral equations: for t > 0 and k = 1, 2, 


X k {f) = X°(t) - N k hF k>e (t) - J (0)D 2 t4(f) - Y k ,(0) 1/2 B 0 , k (F k , e (t)) 

+ f\x k (t - s)) + dF k (s ) + f Ff(t - s)dA(s) - W k (t), (3.42) 

Jo Jo 

Y k (t ) = *fc(f) + N k p k F kie (t) - Dfc(0) 1//2 B 0i fe/(i4;/ )e (t)) + J(0) 1 / 2 (t/ fc (f) - 17(f)) 

- Ax fc (f - s)) + dF k (s) + f F k (t - s)dA(s) + Wf(t) - (3.43) 

Jo Jo 

S(t ) = 5°(f) + D>(0) 1/2 B 0 ,i(74, e (f)) + T 1 (0) 1 / 2 l?o,2(D 2 , e (f)) + J(0) 1 / 2 D(f) + ID(f) + (f), 

(3.44) 

and E , Q , B and D are given as follows: 

E k (t) = A(t) - (X k (t))+, D k (t) = X k ^) + Mt) ~ Mt), (3-45) 
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Q k (t) = (X k (t)) + , B k (t ) = ~(X k (t))~, 

where 

X° k (t) := X k (0)Fl e (t) + (X k (0)) + (Ff(t) - Fg >e (t)), (3.46) 

2 

s°(t) := J2(Xk'(0)F k , e (t) + Z k ,(0)F k {t)F k ,, e (t)) + J(0)F he (t)F 2 , e (t) + J(0)F m (t), (3.47) 

k= 1 

Yg(t) := %{ 0) + X fc (0)F fc , c (*) + (X fc (0))+(F fc (i) - F k , e (t)) - S°(t ), (3.48) 

i/ie processes B o k := {B ok (t) : t > 0}, k = 1,2, are independent standard Brownian bridges, the 
process U is a continuous two-parameter Gaussian process defined above with the processes U\ and 
U 2 defined in (3.40), and the processes W k , Wf and W are defined in (3.35), (3.37) and (3.36), 
and B$ k is independent of U and W k , Wf and W, and the process 4 1 := {4/( t ) : t > 0} defined by 

^(t) := [ f £(t, - si, t - s 2 )dF(s\, s 2 ), (3.49) 

Jo Jo 

is a well-defined continuous process, where, for s\,s 2 > 0, 

<S(si,s 2 ) := £a(si)l(si < s 2 ) + E 2 (s 2 )l(s 2 < si) + {Ei(si) A E 2 (s 2 ))1(si = s 2 ). (3.50) 

We remark that the limit processes X k , k = 1,2, have the same structure as the unique solution 
to an integral convolution equation, as shown in Reed [45], but are also different because they are 
both driven by the same generalized multiparameter Kiefer process K defined in Proposition 3.1. 
These two limiting processes X k , k = 1,2, are correlated because of the correlated service times of 
the parallel tasks of each job, which is captured by the process K , as well as the same arrival limit 
process A. In fact, these two processes K and A as well as the limits associated with the initial 
quantities are the driving stochastic components of all the limit processes in (3.42)-(3.45). 


4 Concluding Remarks and Future Work 

We remark on the main ideas of the proofs for the limit theorems due to space constraint. The 
main difficulty in the study of many-server fork-join networks with NES is the resequencing of 
arrival orders after service completion at each service station. Tasks of distinct jobs must be 
differentiated and tracked in order to describe the waiting buffer dynamics for synchronization. To 
mathematically describe the system dynamics, we develop a new approach using multiparameter 
sequential empirical processes driven by service vectors for parallel tasks of each job, as depicted 
in Figure 2. This approach is used to establish FLLNs and FCLTs for the waiting buffer processes 
for synchronization and the service processes jointly in the fundamental fork-join network where 
all service stations are operating in the many-server heavy-traffic regimes. 

As a prerequisite, we first establish a new FCLT for multiparameter sequential empirical pro¬ 
cesses driven by random vectors (Theorem 2.1). To prove Theorem 2.1, we employ the standard 
approach of establishing convergence of finite-dimensional distributions and tightness [?, 20, 55]. 
The convergence of finite-dimensional distributions follows from the strong convergence result of 
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multiparameter empirical processes in [42], To prove tightness, we present a new decomposition 
property for multiparameter sequential empirical processes, which have a multiparameter mar¬ 
tingale [23, 19], and a second term of finite variation. We apply properties of multiparameter 
martingales [23, 19] and strong approximations of random walks by Brownian motions (see section 
3.5 in [30]) to show the tightness of those two decomposed terms, respectively. This decomposition 
also plays a very important role in proving the tightness of the number of tasks in each waiting 
buffer for synchronization, the number of tasks in each parallel service station and the number 
of synchronized jobs. Specifically, the aforementioned processes can be decomposed into a linear 
combination of three terms: an integral functional of the arrival process and two other terms from 
the decomposition of the multiparameter sequential empirical process. We apply Aldous’ tightness 
criteria (see, e.g., Lemma 3.7 in [28]) and another tightness criteria for processes with proper decom¬ 
positions satisfying certain conditions (Lemma VI.3.32 in [20]) to verify the tightness property of 
the two terms related to the sequential empirical process driven by the service vector, respectively. 

The proofs of the limit theorems in the QD regime can be regarded as generalizations of those 
for G/GI/oo queues in [28]. However, since all the processes, X , Y and S', are represented via 
the multiparameter sequential empirical processes driven by the service vectors, many technical 
challenges must be addressed in the multiparameter setting, for example, using multiparameter L 2 
martingales, and mean-square limits of (integral functionals of) multiparameter processes defined 
on (k > 2). One important advantage of our new approach is that all the diffusion-scale limit 
processes for X, Y and S are all functionals of two independent processes - the arrival limit and 
the multiparameter generalized Kiefer process driven by the service vector (Theorem 2.3). From 
that, the characterization of the joint transient and stationary distributions of these processes is 
made possible (Theorem 2.4). 

The proofs in the QED regime are based on the important observations that the system dy¬ 
namics of G/GI/n queues can be represented via the corresponding G/GI/oo service dynamics 
[45], and that waiting times in the QED regime are 0(1/^/n ) while service times are 0(1). For 
the fork-join network, we represent the dynamics of X, Y and S via that in the corresponding 
infinite-server fork-join network where the entering service times in the model are regarded as the 
“arrival” times for the corresponding infinite-server fork-join network, as shown in Figure 2(b). 
The observation that the entering service times in the parallel stations have a difference of order 
0(1/y/n) is key to prove the joint convergence of the aforementioned processes. On the other 
hand, since we have to simultaneously handle the waiting times of all parallel tasks and work with 
multiparameter sequential empirical processes, we must develop new techniques to prove tightness, 
including establishing new properties for multiparameter L 2 martingales, and identifying a new 
multivariate integral mapping to apply the continuous mapping theorem. 

We believe that a general framework has been developed to study fork-join networks with 
NES in the many-server heavy-traffic regimes (QD and QED). It can be potentially used to study 
performance evaluation, capacity allocation, and control problems in multi-class fork-join networks 
under NES with multi-stage processing. We want to find optimal scheduling and routing policies 
such that delays for synchronization as well as delays for service can be minimized, particularly, 
reducing delays for synchronization to be of a smaller order than service. We also want to find 
optimal staffing policies to stabilize delays for synchronization in addition to delays for service 
when arrival rates are time inhomogeneous. Our methods can be extended to investigate reliability 
of many-server fork-join networks under NES in random environments (e.g., service disruptions). 
Fork-join networks with NES are more likely to suffer from service disruptions due to the structural 
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complexity of parallel and sequential task processing. Component-level unreliability can be much 
more amplified by its large scale. We will extend our approach to investigate the impact of service 
disruptions in one or multiple service stations upon system congestion, particularly, delays for 
synchronization and throughput. 
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